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ABSTRACT 

We use numerical simulations to explore whether direct collapse can lead to the formation of super- 
massive black hole (SMBH) seeds at high redshifts. Using the adaptive mesh refinement code ENZO, 
we follow the evolution of gas within slowly tumbling dark matter (DM) halos of M v ; r ~2x 10 8 M Q 
and i? v i r ~ 1 kpc. For our idealized simulations, we adopt cosmologically motivated DM and baryon 
density profiles and angular momentum distributions. Our principal goal is to understand how the 
collapsing flow overcomes the centrifugal barrier and whether it is subject to fragmentation which 
can potentially lead to star formation, decreasing the seed SMBH mass. We find that the collapse 
proceeds from inside out and leads either to a central runaway or to off-center fragmentation. A 
disk-like configuration is formed inside the centrifugal barrier, growing via accretion. For models 
with a more cuspy DM distribution, the gas collapses more and experiences a bar-like perturbation 
and a central runaway on scales of <, 1 — 10 pc. Wc have followed this inflow down to ~ 10 pc 
(~ 10 AU), where it is estimated to become optically thick. The flow remains isothermal and the 
specific angular momentum, j, is efficiently transferred by gravitational torques in a cascade of nested 
bars. This cascade is triggered by finite perturbations from the large-scale mass distribution and by 
gas self-gravity, and supports a self-similar, disk-like collapse where the axial ratios remain constant. 
The mass accretion rate shows a global minimum on scales of ~ 1 — 10 pc at the time of the central 
runaway. In the collapsing phase, virial supersonic turbulence develops and fragmentation is damped. 
Models with progressively larger initial DM cores evolve similarly, but the timescales become longer. 
In models with more organized initial rotation — when the rotation of spherical shells is constrained 
to be coplanar — a torus forms on scales ~ 20 — 50 pc outside the disk, and appears to be supported by 
turbulent motions driven by accretion from the outside. The overall evolution of the models depends 
on the competition between two timescales, corresponding to the onset of the central runaway and 
of off-center fragmentation. In models with less organized rotation — when the rotation of spherical 
shells is randomized (but the total angular momentum remains unchanged) — the torus is greatly 
weakened, the central accretion timescale is shortened, and off-center fragmentation is suppressed 
— triggering the central runaway even in previously 'stable' models. We estimate the resulting seed 
SMBH masses to lie in the range 2 x 10 4 M Q — 2 x 10 6 Mq, substantially higher than the mass range of 
Population III remnants. Corollaries of this model include a possible correlation between SMBH and 
DM halo masses, and similarity between the SMBH and halo mass functions, at time of formation. 
Subject headings: methods: numerical — galaxies: formation — galaxies: high-redshift — cosmology: 
theory — cosmology: dark ages, reionization, first stars 



1. INTRODUCTION 

Observational evidence points to most large galax- 
ies hosting su permassiye blac k holes (SMBHs) in their 
centers (e.g., iFerraresd 120051 ) and to the possible co- 
evolut ion of galaxies and SMBHs, whether causal 
(e.g.. iGebhardt et all 120001 iFerrarese fc Merrittl l2Cl5(i 
iTremaine et all |2002[) or not (jJahnke fc Macciol 120111 ). 
Recent high-z observations have found bright QSO s 
at z <; 7 (e.g.. iFan etHI [200l iMortlock et all I20TTI ). 
Hence, some SMBHs with M, ^fewx 10 9 M Q had already 
formed before the universe was ~ 700 Myr old, raising 
the issue of how such large black holes (BHs) could have 
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grown so quickly. 

The early formation and growth of SMBHs must be 
understood in a cosmological context. Early SMBHs 
could have formed from remnants of the first genera- 
tion of stars, Population III, whic h subsequently grew by 
gas accretion and mergers (e .g . . lHaiman fc Loebl 2001 



Yoo & Miralda-Escudc 2004; Volontcri & Recs 2006; 
Li et al.ll2007t iPelupessv et all 120071: iTanaka fc Haimanl 
20091 ). Early studies of the initial mass function 



(IMF) of these objects argued in favor of very mas- 
sive M sta .r ~ 100- 1,000 M w objects (e.g. [Abel et all 
120021 IBromm et al.l[2002l IQ'Shea fc Normanl 12001) . but 
more recent wor k indicates a rather normal IMF (e.g., 
iWise et "all 120121) . Even if there were a sufficient sup- 
ply of 100 Mq remnant BHs which grew via gas ac- 
cretion at the Eddington rate with 10% radiation effi- 
ciency, it would h ave taken them at least ~7x 10 8 yr to 
reach ~ 10 9 M Q (|Salpeterlll96l ). This timescale barely 
matches the available time required to make bright QSOs 
at z ~ 6 — 7. 
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Such efficient growth by accretion is unlikely, for two 
reasons. First, Population III stars in the vicinity of 
the seed, and the BH progenitor itself, ionize the sur- 
rounding matter, causing expansion of the HII region 
and likely su ppression of subsequent accre tion onto the 
SMBH (e.g., ISafranek-Shrader etHIIMTl . Second, the 
possible supernova explosion of a Pop III star can evacu- 
ate the ambient gas in the vicinity of the seed. 

This suggests that mergers must play an essential role 
in growing SMBHs from Population III seeds, but this 
is difficult as well. According to the cold dark matter 
(CDM) cosmology, dark matter (DM) halos frequently 
merge in hierarchical structure formation. However, too 
frequent halo mergers would result in the formation of 
small TV-body systems of BH seeds, which would be prone 
to SMBH slingshot cjcction(s) from the DM halo. Given 
the relatively short time available for growing SMBH 
seeds to quasar masses, the long list of delaying processes 
can seriously hurt the Pop III seed scenario. 

An appealing alternative is that early SMBHs formed 
via direct collapse of the ha l o gas at z ~ 10 — 
20 (e.g.. lOhfeHaimanl 12001 iBromm fc Loebl [200l 



Volonteri fe Reesll2005|; IBegelman et al.ll2006UWise et al 

Regan fe Haehnelti 1200 " 



20081: iLevine et al 



2008; 



Begelman fc Shlosmadl2009l:ISchleicher et al.l2010T) . The 
direct collapse paradigm assumes that an SMBH seed, 
M, <; 10 5 Mq, forms from the cold halo gas. This pro- 
cess favors a high density, high inflow rate environment 
in which star formation is suppressed. A massive central 
object — a supermassive star (SMS), forms at the center 
of the collapsing region. This object is powered by a com- 
bination of core nuclear burning and Kelvin - Helmholtz 
contr action (e.g., IBegelman et al.l [20061 120081 : iBegelmanl 
i2010t ). After the stellar core collapses and forms the 
SMBH seed, its convective envelope is powered by the 
accretion onto this seed — this configuration is termed 
a quasistar. Consequently, the seed BH, ~ 100 M e at 
the beginning, can grow to ~ 10 5 Mq in less than a few 
Myrs. 

In this context, constraints on the accretion rate onto 
the SMBH become weaker, and additional problems an- 
ticipated in the Pop III seed scenario do not play a major 
role. In order to make this scenario work, the gas needs 
to inflow rapidly from ~ kpc to ~100AU scales. The 
most probable site for such runaway gas collapse is ex- 
pected in DM halos with virial temperatures ^ 10 4 K. If 
the halo gas can cool only via atomic cooling, the neces- 
sary condition for the collapse to be triggered is for the 
gas temperature to be lower than the virial temperature 
of the host DM halo. 

However, two caveats t o the direct collapse scenari o 
must be addressed (e.g., IBegelman fc Shlosmanl 120091 ). 
First, the angular momentum barrier, in principle, can 
terminate the collapse well before it reaches ~ 100 AU 
scales. Owing to angular momentum conservation, col- 
lapse is halted when the rotational velocity reaches the 
circular velocity. In order to overcome this barrier, there 
should be a physical process to redistribute the gas an- 
gular momentum outward. Second, gas could fragment 
during the collapse, depleting the accretion stream by 
forming gas clumps and ultimately stars. The gas clumps 
could also disturb the accretion pattern. 

In this paper, we study the physical processes that ac- 
company the initial and intermediate stages of gas col- 



lapse in a DM halo. We focus on the dynamical pro- 
cesses: spontaneous and induced symmetry breaking, 
nested gaseous bar formation, and the role of supersonic 
turbulence. Section [2] discusses the relevant processes. 
In Section [3] we describe the numerical setups for a set 
of simulations intended to quantify these processes, and 
in Section U we present our results. We summarize and 
discuss our results in Section [SJ 

2. THEORY 
2.1. Angular Momentum Transfer 

It is generally understood that local viscous trans- 
port mechanisms are inefficient on galactic scales. 
Even on scales of a few parsecs, their characteristic 
timcscalcs are prohibitively long and angular momen- 
tum transport therefore requires non-local mechanisms, 
e.g., generated by large-scale magnetic fields. On larger 
spatial scales, the angular momentum redistribution 
most likely depends on long-range gravitational torques 
(e.g iLvnden-Bell fc KalnaisH972tlSrilosman et~aT1ll989l 
I1990T ). When angular momentum is conserved, gas 
collapse is quickly stopped by the centrifugal barrier. 
Prior to collapse within DM halos, the gas can ex- 
hibit the same angular momentum distribution as the 
DM. The typical spin acquired by a halo during max- 
imal expansion is given by a dimensionlcss parameter 
A = J/v2-Mvi r u c i? v i r , where J is the halo angular mo- 
mentum, and Af v ; r and v c are its virial ma ss and the cir- 
cular velocity at the virial r a dius itW (e.g., lPeeblesiri969l : 
IBarnes fc Efstathioul fl987l : Bullock et all 120011) . The 
mean spin of DM halos is A ~ 0.05. If the gravitational 
potential of the gas dominates, this A will bring the gas 
to a centrifugal barrier when it has collapsed by a factor 
of ~ 100. When the DM potential dominate s, the col- 
lapse will stop after a decade in radius (e.g., I Mo et all 
H998tlShlosmar]l2013l and refs. therein). 

Inflow beyond the centrifugal barrier requires substan- 
tial and continuous angular momentum loss by the gas. 
What are the options? If the gas disk becomes non- 
axisymmetric, gravitational torques between the inner 
and outer gas, as well as between the gas and DM, can 
drain angular momentum away from the collapsing gas. 
The lowest non-axisymmetric modes, m = 1 — 3, pre- 
vail in both the gas and the DM, because their rise 
times are shortest. The m = 1 mode requires pertur- 
bation of the center of mass of the gaseous disk — this 
is not always possible. The next fastest growing mode 
is 771 = 2 — the bar-like mode. Typically it will have 
the highest amplitude, in the absence of 777 = 1. For 
the inflow to extend to smaller spatial scales, the gas 
disk must be self-gravitating. Under these conditions, a 
cascade of bars can be maintained over a substantial dy- 
namic range in radius. Contraction of the gas via such 
an avalanche can be related to s elf-organized criticality 
(|Lichtenberg fc Liebermanl |1995l ) . An analogy can be 
made between the chaos driven when the bar exceeds a 
certain strength and a s and pile whose sl ope is increased 
until the sand slides off (Shlosman 2005). 

The bar cascade can be triggered in a number of ways. 
First, if the disk becomes self-gravitating and cold, its 
axial symmetry is known to be broken spontaneously 
when the ratio of bulk kinetic energy to absolute poten- 
tial energy, T/|W|, is larger than a certain critical value. 
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The gas bar exerts torque on the gas disk and transfers 
angu l ar momentum outward (e.g., jShlosman et al.lll989l 
119901 : lEnglmaier fe Shlosmanll2004l) . This process allows 
the collapsing gas to cross the centrifugal barrier and to 
continue the collapse to smaller scales, if angular mo- 
mentum is continuously extracted. The threshold value 
for gaseous bar formation depends also on the shape 
of the gas distribution an d is a = 0.5fT/\W\ > 0.34 
(|Christodoulou et al.l 119951) . where / = 1 for disks and 
/ = § for spheres. This instability is spontaneous and 
does not require any trigger — it will grow exponentially 
from any infinitesimal perturbation. 

Several additional effects can drive a cascade through 
finite amplitude perturbations. First, the shape 
and dynamics of a self-gravitating gas disk will re- 
spond to the potential of the DM halo in which 
it is embe dded. DM halos are known to be tri- 
atrial (e.g., lAllgood et all 12 006; Bc rentzen et aLl 120061: 
iBerentzen fc Shlosmanl 120061: iHeller et al.l I2007J) . and 
therefore exert gravitational torques on the disk, which 
will respond with low Fourier modes. Second, finite 
bar perturbations can also be driven by tidal effects 
of massive DM and baryon clumps, i.e., halo substruc- 
ture (e.g., iRomano-Diaz et al.l 120081) and from mergers. 
In addition, asymmetry in the background gravitational 
potential can be generated by the spontaneous break 
of axial symmetry in the collapsing flow, DM, gas or 
both. It is well-known that non-radial perturbations 
can grow i n supersonic accretion flows onto compact ob- 
jects (e.g., iHunteii fl96l iGarlickl H97§ iMoncriefl [l980t 
iGoldreich et al.l 119961: lLai fe Goldreicbi 12000J), without 
any assistance from self-gravity. Gas collapse inside an 
axisymmetric DM halo is similar — in Section |4~T1 of the 
present work, we show the appearance of non-radial per- 
turbations in the equatorial plane of the disk, although 
we start from idealized conditions of an isolated halo and 
embedded baryons which possess axial symmetry. In a 
comprehensive cosmological simulation, all these factors 
would provide finite amplitude perturbations that do not 
require the a-parameter to cross the threshold. 

At the spatial scale on which the gas disk forms, the 
potential contribution from the gas is already compara- 
ble with that of the DM, and the gas cooling time is 
shorter than the dynamical time. These conditions im- 
ply that the collapsing gas becomes self-gravitating and 
nearly isothermal, and the density profile evolves toward 
a singular isothermal spherql p ex R~ 2 . Before the gas 
reaches the centrifugal barrier, we can assume it flows 
in with the free-fall speed, whose dependence on radius 
is weaker than logarithmic. Away from the disk and 
the centrifugal barrier, the gravitational potential is still 
dominated by the DM, so the accretion rate can be esti- 
mated at 
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where v 30 = vg/SOkms -1 for a 2 x 10 8 M© DM halo, 
and we used the universal baryon fraction for the gas-to- 
total mass ratio. At these radii, only the gas participates 
in the collapse, as the DM is in equilibrium due to the 
initial conditions. 

4 Throughout this paper, we use _R for spherical radii, and r for 
cylindrical radii. 



On the other hand, when the gas dominates the po- 
tential, the mass accretion rate is of order ~ a v c^/G ex 
T 3 / 2 /G, where c s is the gas sound speed and a v is 
the viscosity parameter in the disk (e.g. iShul 119771 : 
iShlosman fe Begelmanl 119891 ) . As we shall see in Sec- 
tion 14- 1 1 this situation in our simulations occurs in- 
side the centrifugal barrier. For local viscosity in the 
disk, a v < 1, but for non-local viscosities, e.g., mag- 
netic or gravitational torques, a v ^ 1, especially in 
the self-gravitating c ase where the effective a v ^> 1 
(|Shlosman et al.lfl990f) . In this latter case, the inflow rate 
is still given by Equation Q] where vs is the virial veloc- 
ity of the gravitational potential which is now dominated 
by the gas, i.e., M gas /M to t = 1. The inflow velocity will 
depend on the compactness of the gas distribution and 
can substantially exceed the virial speed of the DM halo. 

A high accretion rate is crucial for the formation of 
a SMS or di sk that can dev elop into an SMBH seed. 
For example. iBegelmanl (|2010f ) pointed out that an infall 
rate exceeding ~ IMq yr _1 is necessary in order to form 
a 10 6 M Q SMS, which has a nuclear burning timescale 
of only a few million years. The initial SMBH seed 
of 100 — 1,000M Q , formed by core collapse of such a 
star, then grows at a highly super-Eddington rate inside 
the remaining convective envelope (the 'quasistar' phase: 
iBegelman et al.l (|2008f )). reaching ~ 10 5 Mq or more in 
less than a few Myr. 

2.2. Fragmentation of Accretion Flows 

Gas fragmentation can terminate the collapse by con- 
suming the gas supply. In addition, the clumps formed 
during fragmentation can excite odd-mode perturbations 
in the gas disk, damping the bar instability which plays 
the key role in angular momentum redistribution. In 
order to continue collapse to very small spatial scales, 
fragmentation should be suppressed. 

There are several ways to suppress gas fragmenta- 
tion. Supersonic turbulence can suppress fragmenta- 
tion via shocks that sweep the density fluctuations in 
a crossing time (e.g.. lPadoanlll995l ; IPadoan fe Nordlundl 
120021 : iKrumholz fe McKed l2005f ) . Studies of supersonic 
turbulence find that the probability distribution func- 
tion (PDF) of the density fluctuations is lognormal (i.e., 
Gaus sian in the log) under a wide range of conditions 
(e.g., lyazquez-Sema deni 19 941; IPadoa n 1995; Scalp et all 
199a lOstriker et all Il999t IPadoan fc Nordlundl 120021: 



Krumholz fc McKec 2005): 
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where the distribution mean In 

(po being mean density), and its dispersion is o* ~ 

[ln(l + 3.M 2 / 4 )] (M is the flow Mach number). This dis- 
tribution is the result of x being a random variable, which 
is itself a product of independent random variables. Note 
that the lognormal density PDF has been observed in 
simulations of turbulent motions on the scales of giant 
molecular clouds, where the background (unperturbed) 
density is uniform. Here, we shall test the development 
of turbulence in collapsing flows with steep preexisting 
density gradients. 
The critical overdensity required for fragmentation is 
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x cr it = (Ajo/A s ) 2 , where Ajo is the Jeans length at the 
mean density and A s is the turbulent length scale at 
which the velocity d i spersi o n reaches the sound speed 
(iPadoan fc Nordlun d 2002; iKrumholz fc McKed [20051 : 
IBegelman fc Shlosm an 2009). In a supersonic turbulent 
medium, overdensities higher than x cr ;t are unstable to 
fragmentation. Estimating the mass fraction in Jeans- 
unstable fragments that collapse on a timescale shorter 
than the free-fall timescale as / = f°° x p ^' dx, we find 
that / < 0.02 for A4 > 3. If the collapsing gas maintains 
supersonic turbulence, the amount of fragmenting gas is 
negligible. 

In our models, we have assumed optically thin flows. 
The spherically-symmetric, isothermal density distribu- 
tion of baryons modeled in this paper becomes opaque at 
~7x 10~ 3 A pc, where X is the ionization fraction. For 
temperatures encountered in this simulation, X<Cl and 
the electron scattering optical depth is negligible, as is 
the free-free absorption. For primordial chemical compo- 
sition and the range of temperatures ~ 3, 000—10, 000 K, 
the opacity is expected to be dominated by bound-bound 
transitions in hydrogen, and especially by the opacity for 
Lyman a photons. 

Radiation fields can suppress gas fragmentation by dis- 
sociating molecular gas or preventing its formation in 
the first place. If H2 is absent, the g as temperature 
will be maintaine d at > 8,000K (e.g., lOmukail l200ll: 
IShang et al.l I2010D . The prime agent of H 2 dissocia- 
tion can be an externally-produced Lyman- Werner (~ 
11 eV— 13 eV ) continuum, e . g., fro m neighboring Pop III 
stars, (e.g.. iDiikstra et aLl 120081 ). H 2 formation could 
be inhibited by the high temperatures likely to obtain 
in the collapsing gas. The cooling efficiency of Lya pro- 
duced in the accretion flow at T ~ 10 4 K is limited by 
the high optical depths. For a sufficiently high column 
density of neutral hydrogen, collisional de-excitation can 
decrease the escape fraction of Lya, preventing the gas 
from cooling below ~ 8, 000K and suppressing H2 forma- 
tion thermally (e.g.JSpaans fc Silk 2006; Schleicher et al.l 
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One might speculate whether Lya photons generated 
within the accretion flow could be upscattered into the 
Lyman- Werner continuum before escaping, thus produc- 
ing a dissociating flux self-consistently. For an isother- 
mal density profile, we estimate the required Lya optical 
depth (measured in the Doppler core) of r Q ;> 10 14 (e.g., 
Harrington! 1 19731; iNeufeldl 119901 ). which can be attained 
at radii <, 10~ 3 pc for the parameters of our models. 
Other processes, however, are likely to allow the energy 
in these upscattered Lya photons to leak away before 
they can escape. For one thing, these photons can be 
absorbed by ev en small amou nts of dust. Extrapolat- 
ing Figure 18 of INeufeldl (J1990D to our column densities 
(which can be inferred from the next section), we find 
that a metallicity as low as 10 -8 solar will do the job. 
An additional process that can destroy the Lya photons 
is the resonant pumping of H2, followed by fluorescent 
decay to vibrational levels. 

We note that fragmentation might be avoidable in the 
prese nce of supersonic turbu l ence, even if H2 is present 
(e.g., IBegelman fc Shlosman! [20091) . This happens be- 
cause the fraction of fragmenting gas decreases with M . 

In our simulations, we assume that H2 is destroyed and 



do not analyze its contribution to the gas cooling. We 
also neglect magnetic fields and their effects on the turbu- 
lent flow, and model gravitational collapse within the DM 
halo hydrodynamically. This is justified because i?-ficlds 
are expected to be weak at high z, and because of the 
high temperature of the flow, ^ 1,000K. At these high 
temperatures, the value of the critical B-field strength 
needed to limit the co mpression in supers onic isothermal 
shocks is also higher ([Padoan et al.112007}) . 

3. NUMERICAL TECHNIQUE 
3.1. The numerical code ENZO 

In order to test the theoretical arguments made in Sec- 
tion [21 we use an Eulerian adaptive mesh refinement 
(AMR) code ENZO-2.1, which has been tested exten- 
sively and is publicly a vailable ([Bryan fc Normanl 119971 ; 
iNorman fc Brvan|[l999T ) . ENZO uses a particle-mesh N- 
body method to calculate the gravitational dynamics in- 
cluding collisionless DM particles, and a second-order 
piecewise parabolic method (PPM, iBrvan et al.1 119951 ) 
to solve hydrodynamics. The structured AMR used in 
ENZO places no fundamental restrictions on the number 
of rectangular grids used to cover some region of space 
at a given leve l of refinement, or on th e number of levels 
of refinement (jBerger fc Colellal 119891 ). A region of the 
simulation grid is refined by a factor of 2 in length scale if 
the gas density is greater than poN 1 , where po is the min- 
imum density above which refinement occurs, N = 2 is 
the refinement factor and I is the AMR refinement level. 
The iTruelove et ail (|1997|) requirement for resolution of 
the Jeans length is verified. 

ENZO follows the non-equilibrium evolu tion of six 
species: H, H + , He He + , He ++ , and e~ (jAbel et al.l 
119971 ; lAnninos et al.|[l997l) in a gas with primordial com- 
position. It calculates radiative heating and cooling fol- 
lowing atomic line excitation, recombination, collisional 
excitation and free- free transitions. Radiative losses from 
atomic cooling are computed in the optically-thin limit. 
As we discussed in Section [21 we exclude the chemistry 
and cooling related to H2 in this paper. 

3.2. Simulation Setups 

For a gas with primordial composition and no H2, the 
earliest collapse can occur in DM halos whose virial tem- 
peratures exceed ~ 10 4 K. For this reason we set up 
a spherical DM halo with M vir ~2x 10 8 /i _1 M© at 
z = 15. According to the top-hat model, such a halo 
will have a virial radius of i? v ir ~ 945/i -1 pc and a virial 
temperature T vir ~ 32, 000 K, according to the WMAP5 
cosmology (£U= 0.2 79, fi b = 0.0445, and h = 0.701) 
(jKomatsu et al.ll2009f ) . We simulate this DM halo within 
a (6kpc) 3 computational domain. 

For the DM halo, we assume an isothermal sphere. 
This halo mod el is similar to the univer sal DM halo 
model (NFW, iNavarro et all 119961 I1997D at the spa- 
tial scales of interest. The isothermal model also pro- 
vides a simple analytical form for the velocity distribu- 
tion function which allows us to generate a live isotropic 
DM particle distribution that maintains a stable equi- 
librium with 10 6 particles. We introduce a DM core 
radius, Rdm, within which the DM density is approxi- 
mate l y constant (e.g.. [El-Zant et al J l200.lt iTonini et al.l 
l2006t iRonTano-Diaz et al.1 120081; iPrimackl I2009D . There- 
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TABLE 1 
DM Core Radius and Central Collapse Time 



Models # DM Core Radius _R dm (pc) Collapse Time (Myr) 



A 
B 
C 
D 

E 



0.10 
0.40 
0.75 
1.50 
6.00 



2.1 

1.7 

8.7 
no collapse 
no collapse 



Dmod 



1.50 



1.5 



Note. — The collapse time given in the third column is time 
between the start of the simulation and the start of the central 
runaway collapse. 

fore, i?dm plays the role of the King radius for a nonsin- 
gular isothermal sphere. We vary i?dm to see the effect of 
DM halo structure on the development of gas collapse. In 
particular, we implement a number of simulations, mod- 
els A - E, with different DM density core sizes, as given 
in Table [TJ 

The DM halo and the gas have been laid down at virial 
and pressure (for the gas) equilibrium, so no transient ad- 
justment occurs. The gas halo is taken to be an isother- 
mal sphere as well, with an R COIC = 100 pc constant 
density core at the start of the simulations. The total 
gas mass in the halo is estimated from the cosmological 
baryon fraction, J7b/^ m ~ 0.16. The angular momentum 
of the halo gas is computed from J = Av^Mviri'c-Rvir, 
defined in Secti on El In large-scale iV-body cosmologi- 
cal simulations, iBullock et al.l (|2001|) has found that the 
specific angular momentum in DM halos roughly follows 
j(R) ~ R 11 , close to a flat rotation curve. We therefore 
impose a constant rotation velocity on the gas outside the 
core, assuming that the rotational axis runs parallel to 
the z-axis. For R <, i? C ore, we assu me solid body rotation 
for the gas, i.e., i(R) ~ R 2 (e.g., iSamland k, Gerhardl 
l2003t iHelier et alll200l . 

The outer boundary of the model is treated in the fol- 
lowing way. We smooth the sharp density cutoff at the 
virial radius by adopting the gas density profile ~ R~ 3 
beyond i? v ir- As the gas is not in exact pressure equi- 
librium at this point, this leads to its slight expansion 
which does not affect the simulations. 

We also ran a number of models with modified initial 
conditions. One such representative model is discussed 
as model Dmod, i.e., modified model D, in Section l4~4l 
In this model we randomized the orientation of the ini- 
tial angular momentum in spherical shells, keeping the 
mean-square angular momentum unchanged. For this to 
happen, we have slightly decreased j of the inner shells 
and compensated the outer shells by a slight increase in j. 
This was done in order to mimic cos mological initial con- 
dition s, based on the simulations by iRomano-Diaz et al.l 
(|2009T) . in particular on their Figure 19. 

To verify that the DM halo is indeed formed in equi- 
librium, we tested DM-only models and confirmed their 
stability, especially in the central region. Additional tests 
have been run for model B with an isothermal equation 
of state for the gas; these show very similar evolution. 
We also tested a model with an adiabatic equation of 
state — the gravitational collapse did not proceed far in 
this model, as expected. 



4. RESULTS 

We would like to emphasize several important points 
about the gravitational collapse modeled here, before 
we show the results of numerical simulations. First, 
the evolution described here truly proceeds from inside 
out. Development of the central mass accumulation and, 
therefore, of the prospective seed SMBH, happens on 
timescalcs < 10 Myr, much shorter than the free-fall 
timescale for the DM halo, ~ 80 Myr. Hence, the dy- 
namics of the outer parts on scales of ~ 0.1 — lkpc, 
although interesting in itself, appears to be irrelevant 
for the central regions of <, 10 — 50 pc, which dominate 
the formation of the SMBH seed. Next, the gravita- 
tional collapse proceeds in two stages. The first stage 
involves infall, braking and the formation of a gaseous 
disk-like configuration inside the centrifugal barrier. All 
the models exhibit this phase. Models A - C show the 
second stage of gravitational collapse in which the inner 
part of the disk, r <, 1 — 5pc, develops a runaway which 
proceeds to the point where the numerical evolution has 
been terminated, at ~ 10 -4 pc ~ 20 AU. Our task will be 
to explain this evolution. We shall emphasize models B 
and D as representative of 2-stage and 1-stage collapse, 
respectively. Model Dmod is discussed separately. 

4.1. Loss of axial symmetry and central runaway 

As the gas is not fully supported rotationally, it goes 
into nearly free-fall collapse, which develops from inside 
out because of the shorter dynamical timescales at small 
r. There is little difference among models at this stage, as 
they differ only in the value of the DM core radius i?dm , 
and, therefore, in the depth of the DM potential well. 
Figure [T] shows face-on density- weighted projections of 
the densitjo of the gas disk at t = 4.6 Myr in models 
B and D, with i?dm ~ 0.4 pc and 1.5 pc, respectively. 
Unless mentioned, other models behave similarly. 

During the first stage of collapse, the baryon angu- 
lar momentum flows inward. The inner gas is the first 
to reach the centrifugal barrier, and develops a disk-like 
configuration at r ~ 1 pc. The disk grows quickly in ra- 
dius and mass, and by t ~ 4 Myr, the snapshots display 
well-developed gas disks with radii ~ 10 pc. The disk 
boundaries, both in r and in z, are delineated by stand- 
ing shocks, discussed later. The surface densities of all 
models are well-approximated by a power-law S ~ r~ n 
with n ~ 1.3. Figure [2] displays the evolution of the sur- 
face density for model B, as well as comparisons of the 
surface densities at t = 4.7 Myr among models B, D and 
E. The disks are truncated at r ~ 10 pc, as clearly seen 
in this figure. The outer surface density profiles of all 
disks, in models A to E, are very similar, but in the cen- 
tral ~ 1 — 2 pc the density profiles differ: the density in 
model B within the central 0.1 pc is more than an order 
of magnitude higher than the density in model D. 

The DM density profiles show little evolution even in 
the central regions. In models B, D and E (Figure [2]), 
and also in all other models, we observe the trend that 
more cuspy DM halos develop higher surface and volume 
density gas disks. The disk in D has a surface density 
about a factor of ~ 10 times higher than E. The disk vol- 
ume density also appears higher in more cuspy DM halos 

5 Ei(piWi)/Si(Wi), where the weight function W\ is p 
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Fig. 1 . — Density-weighted projection of density (see text) of (a) model B, and (b) model D, with thickness Az = 2 pc about the equatorial 
plane, at t = 4.6 Myr. These face-on views show gas disks that have formed from the collapsing gas, immediately after the gas disk in 
model B has undergone runaway collapse at its center. Model B shows the formation of a gas bar that redistributes angular momentum and 
drives a strong mass inflow, whereas the gas disk in Run D maintains a steady density profile. Both gas disks are intermittently turbulent. 
Snapshot resolutions are 0.01 pc (left) and 0.1 pc (right). Box size is 16 pc. 













~ 


a) 










- 






"* "^Trr"-—- 






- 












V 




t= 4.0 Myr 


- - - 


t = 4.6 Myr 








■\ 


— 


t= 4.7 Myr 








•\ 
















A ' 


1 


1 




- 


b) 










- 










^.—^ 


■^-^^^ 


- 











Model B 












" ~ " 


Model D 
Model E 










V 















l°gl( 



-1 
[pc] 



Fig. 2. — Surface density of the growing disk, averaged over 
annuli, during the central runaway, (a) model B (solid blue) at 
t ~ 4.0 - 4.7 Myr. (b) models B (solid blue), D (dashed red) and 
E (dotted black) are compared at t = 4.7 Myr. 



(Figure [3]). This difference is clearly observed within the 
DM core radius i?dm of less cuspy halos. 

Over the next time period of a fcwxl0 5 yr, the disk 
surface density grows as a result of accretion. The sec- 
ond stage of the collapse is reached only in models A 
- C, and is characterized by both baryonic angular mo- 
mentum outflow and continuing mass influx. On larger 
scales, just outside the centrifugal barrier, j increases, 
moving the barrier outward. Inside the centrifugal bar- 
rier, the disk develops a global instability, dominated by 
the Fourier component to = 2 — a bar-like mode, which 
transfers j to the outer gas and to the DM. For exam- 
ple, model B exhibits a well defined gaseous bar in the 
central ^ 1 — 2 pc of the disk at ~ 4.7 Myr. The appear- 
ance of this bar-like mode is characteristic also of smaller 
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Fig. 3. — Evolution of gas volume density profiles averaged over 
spherical shells, for models B and D. Both simulations start from 
identical initial conditions. Owing to the deeper DM potential, the 
gas disk in model B has a steeper mass distribution in the center 
(see the profiles at t = 4.0 Myr). At t = 4.7 Myr, model B shows 
runaway collapse (strong mass inflow) driven by gaseous bars, while 
the model D disk shows marginal evolution from t = 4.0 Myr to 
t = 4.7 Myr. The density profile of the runaway collapsing gas 
structure is close to that of a singular isothermal sphere (dashed 
line). 

spatial scales of this model (Fig. 2]), as we shall analyze 
below. The timescale for the onset of this stage of the 
central runaway is given in Table [TJ and is increasing 
along the sequence A— > C, i.e., toward less cuspy halos. 
On the contrary, models D and E display a weak central 
depression which has the appearance of a ring-like struc- 
ture around the center and no to = 2 axial symmetry- 
breaking. At a later stage, these models exhibit off-center 
fragmentation. 

Figure [TJ, confirms the crucial role of gravitational 
torques in shaping the inner mass distribution via an- 
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Fig. 4. — Face-on density-weighted projection of density showing the gaseous bar cascade in model B on different spatial scales. Upper 
left: f6pc box at t = 4.6 Myr; Upper right: lpc; Lower left: 0.1 pc; Lower right: 0.01 pc at t = 4.7 Myr, the time of the runaway collapse 
at the center. The slice thickness Az is (a) 2pc, (b) 0.5 pc, (c) O.lpc and (d) 0.01 pc. 



gular momentum transfer, which requires a substantial 
departure from axial symmetry. This symmetry breaking 
gives rise to the lowest modes, m = 1,2, and occasionally 
m = 3. Only m = 1 exhibits a mild displacement of the 
center of mass of the disk, while m = 2 has the appear- 
ance of a strong gaseous bar which drives spiral structure, 
and m = 3 shows up as tri-armed spirals. We follow the 
runaway collapse to a spatial scale of ~ 10~ 4 pc (20 AU), 
and stop the calculation there because the flow is ex- 
pected to become opaque and enter a dynamical regime 
where radiation pressure must be taken into account (see 
Section |2~2|) . 

Owing to the AMR nature of ENZO, the small scale 
structure in Figures [2] and [3] is only resolved when a 
given region exceeds the density threshold for a new re- 
finement. The density profile in model B is resolved 
to ~ 20 AU at the end of the simulation, while model 
D gas only reaches a resolution of ~ 0.1 pc. Resolv- 
ing the ~ 20 AU scale means that the gas inflow has 
reached this scale by overcoming the angular momentum 
barrier. It demonstrates that the gravitational torques 
resulting from axial symmetry-breaking at ~ 1 pc in 
model B trigger continuous gas inflow on progressively 
smaller scales. Finally, Figure [3] demonstrates that the 



timescalc of this inflow is very short, from t = 4.0 Myr to 
t = 4.7 Myr, a truly runaway collapse with a dynamical 
time of Ss 10 6 yrs. 

For the idealized case of spherical accretion within an 
external gravitational potential given by p ex R~ 2 , the 
gas mass flux is given by M ~ R 2 pvn with ur weakly de- 
pendent on R. Figure [3] confirms that the density profile 
of the collapsing gas in model B tends to the isothermal 
density profile, p ex R~ 2 , for R ^ 10 pc. It decreases 
sharply outside this radius due to the shock at the cen- 
trifugal barrier, then continues with the same slope at 
larger radii. (Note that all models start with a flat gas 
density core of 100 pc.) Figure [S^, shows the radial pro- 
file of the inflow velocity measured in the disk plane. We 
can divide the r-rangc into the collapsing gas in the re- 
gion outside 10 pc, the disk plateau at r ~ 1 — 10 pc, 
and the inner region of runaway collapse — all at the 
time of the central runaway, t ~ 4.7 Myr. The stronger 
than expected increase of Mach number with decreas- 
ing radius has its origin in a combination of temperature 
decrease (due to an increased gas density) and the steep- 
ening increase of v-r, toward the center, as the result of 
the runaway collapse of a finite gas mass. 

Thus the outermost region is dominated by free-fall 
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Fig. 5. — (a) Evolution of the radial velocity profile as a function 
of r in the disk plane for model B, normalized by the sound speed 
(i.e., in Mach numbers), (b ) Evolution of the azimuth-averaged 
tangential velocity (in Mach numbers) as a function of r in the 
disk plane for model B. 
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Fig. 6. — Evolution of the mass accretion rate in model B at 
various times, from before the second stage of collapse at t = 3 Myr, 
through the runaway at t ~ 4.7 Myr. 

kinematics because the angular momentum is low. The 
density profile is maintained as p ~ R~ 2 if one neglects 
the weak variation of vr. The intermediate disky re- 
gion has ur ~ const, and the same —2 slope of log p. 
The innermost collapse exhibits a rapidly increasing and 
steepening infall velocity (Fig. |SJi,b). The accretion flow 
is mostly rotationally supported in the disky region of 
1 — 10 pc, as well as in the central region, where the 
self-gravitating collapse proceeds with a non-negligible 
angular momentum. The time evolution of vr in the 
central region reflects the self-gravitating collapse devel- 
oping there, and a dynamical decoupling of the gas from 
the DM background potential. 
The measured mass accretion rate, M(R), appears to 
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Fig. 7. — Density profiles along the axes noted at the runaway 
collapse time t ~ 4.7 Myr for model B. The profile along the z-axis 
has been measured at x = y = 0. Note the positions of the radial 
shock at r~ 10 pc and the surface shock at z ~ 0.5 pc. 

be approximately flat, ~ 0.1 — 0.2M Q yr _1 , exterior to 
the radial shock, in agreement with Eq. Q] (Fig. [6]) . In the 
disk region, M(R) decreases abruptly to ~ 0.01 M Q yr _1 
at t ~ 3 Myr, and drops by an additional order of mag- 
nitude during the second stage of the collapse, when the 
material is drained by the central runaway. At this latter 
time, t = 4.7 Myr, the central runaway is accompanied 
by peak values of M(R) ~ 2M Q yr _1 . The decrease in 
the accretion rate in the range of r ~ 1 — 5 pc is a reflec- 
tion of mass conservation, as the bar-like mode channels 
gas inward at a higher rate than it can be resupplied by 
disk accretion, and thus drains this region. This mini- 
mum in M(R) is found in the region where the turbulence 
induced by the radial shock decays. In Section 15.11 we 
use this radius in model B to estimate the baryon mass 
that participates in the second stage of the collapse, and 
therefore, the mass of the seed SMBH. 

The high mass accretion rate is an important ingredi- 
ent o f early SMBH formation in the direct collapse model 
(e.g iBegelman et al|2006l . 120081: iBegelman fc Shlosmanl 
l2009t lBegelmanll2010|) . The central runaway region ex- 
ceeds the mass accretion rate given in Eq. Q] by about 
an order of magnitude. The reason for this lies in the 
fact that during the second stage of the collapse, the 
inflow velocities are determined by the compactness of 
the gas distribution (Section I2.1[) and not by the DM 
halo virial velocity. Under these conditions, the ratio 



M gas /M t , 



1 in Eq. [1] and the free-fall velocity, vg, 



becomes a rapidly growing function of time. As a result, 
the characteristic r-profile of radial infall velocity during 
this stage is characteristic of the 'avalanche' behavior of 
a dynamically decoupled gas. This radial velocity is ex- 
pected to grow sharply with time. The timescale of this 
process depends only on the mass involved in the col- 
lapse, and on the ability of the gas to cool ahead of the 
free-fall time. If the latter condition is fulfilled, the gas 
will collapse to an infinite density in a finite time. 

We note that the gas joining the disk experiences 
strong radial and vertical (z-axis) shocks, as is seen in 
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Figure [3^ (and Figs. [7] and [T3")) . But there is no asso- 
ciated shock in the azimuthal motion. The face-on and 
edge-on velocity fields on scales of ~ 10 pc are shown in 
Figure [Jjjl While the face-on view is dominated by rota- 
tion, the edge-on view exhibits turbulent motions domi- 
nated by vortices, as we discuss in the next section. 

We next examine the geometry of the collapsing gas 
during the runaway stage: is it best approximated as one- 
dimensional collapse with a spherical symmetry, where 
the angular momentum is completely unimportant, or 
does it exhibits a cylindrical or disk-like geometry? To 
answer this, we compiled the density profiles of the col- 
lapsing gas in model B at t ~ 4.7 Myr, along two direc- 
tions in the equatorial plane, ±x and ±j/, and along the 
±z-axis (Fig. [7]). We verify that the collapse outside the 
centrifugal barrier proceeds in a spherical fashion, which 
is understandable because the angular momentum in this 
region is not important. At t ~ 4.7 Myr, the radial shock 
in the equatorial plane is positioned at r ~ 10 pc (see 
also Figure [5^,). 

Inside the radial shock marking the centrifugal barrier, 
i.e., between r ~ lpc and 10 pc, the equatorial density 
increases substantially while the density along the z— axis 
is almost constant. At ~ lpc, the density along the z- 
axis is ~ 10~ 3 times the density measured along the x- 
or y-axis at the same radial distance. The reason for 
this abrupt change in density ratio lies in the relative 
positions of the radial shock and the surface shock (i.e., 
the shock in the vertical velocity, at roughly constant 
z) in the disk. While the centrifugal barrier stops the 
cold inflow toward the rotation axis at r ~ 10 pc, the 
position of the disk surface shock, i.e., the disk thickness, 
is determined by the postshock temperature which never 
exceeds ~ 10 4 K, as seen in Fig. [9] 

A strong surface shock is maintained by the infall along 
the z-axis at ~ 0.5 pc (Fig. [7]). Both radial and surface 
shocks are slowly moving outward with time. Most in- 
terestingly, at the time of the central runaway, the sur- 
face shock collapses toward the equatorial plane around 
x = y = 0. The 'peanut' shape of the surface shock at 
this time can be inferred from Fig. 1131 Figure [7] reveals 
that within the central region the second stage of the col- 
lapse proceeds in the self-similar fashion, preserving the 
disk-like configuration. We have tested this by plotting 
these distributions at various times. But even as a single 
snapshot, Figure [7] nevertheless reflects the evolutionary 
trend because different r are associated with different 
dynamical timescales. The fact that p(z)/p(r) ~ const. 
for different r means that this ratio also does not evolve 
in time — a clear sign of self-similarity during the cen- 
tral runaway. Hence the gas configuration remains dy- 
namically cold, sustaining the bar cascade that efficiently 
transfers angular momentum outward (Fig. |4j. 

We now return to Figure |4l which displays the non- 
axisymmctric bar-like perturbations on a number of spa- 
tial scales. We first observe this instability develop- 
ing on scales ~ 1 — 2pc, and then propagating inward. 
The dynamics of gase o us ba rs has been analyzed by 
lEnglmaier fc Shlosmanl (|2004h . who focused on the de- 
coupling of a gaseous bar from large-scale stellar bars. 
This decoupling is associated with the rapidly increasing 
pattern speed of the gaseous bar and its radial contrac- 
tion. 

We measured the pattern speed of the gaseous bar in 
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Fig. 8. — Evolution of the Fourier amplitudes of modes with 
m = 1, 2, and 3, normalized to the amplitude for m = 0, for 
r < 1 pc and Az = 0.25 pc about the equatorial plane. 

model B on scales of 1 pc and and 0.01 pc. The resulting 
values of the pattern speeds, fibar, confirm that the m = 
2 mode tumbles with substantially different speeds. On 
~ lpc scale, fibar ~4x 10 _5 yr _1 , while on ~ 0.01 pc 
scale, f^bar ~ 5 x 10 _3 yr _1 , about 2 orders of magnitude 
faster. 

While these pattern speeds are well-defined, there is 
no clear 'material' separation between the spatial scales, 
i.e., no separation between bars. This continuity of the 
flow properties is associated with continuity of pattern 
speeds -- smaller scales correspond to larger pattern 
speeds, fibar(0 ~ r~ a , with a ~ 1. The parameter a 
is determined by the mass distribution M(r) inside the 
collapsing gas at the onset of the central runaway. The 
gas volume density scales as p ~ r~ 2 in the central region 
(Fig. [7]) and dominates over the DM density distribution 
there, causing the gas to decouple dynamically from the 
DM background in the region of the central runaway. 
The pattern speed is given by fibar(?") ~ [M(r)/r 3 ] 1 / 2 , 
where M{r) ~ r, which explains the instantaneous value 
of a above. 

To verify that the bar-like (m = 2) mode dominates 
over other modes, e.g., m = 1 and 3, we have Fourier 
analyzed the gas response to the asymmetric potential. 
The dominant mode at early times is the m = 1 mode — 
the disk center of mass is initially perturbed by the asym- 
metry developed in the overall mass distribution (Fig. [5]) . 
This mode, however, decays quickly to a negligible am- 
plitude. An additional odd mode, m = 3, is also present, 
although at lower amplitude — it decays similarly. On 
the other hand, the even mode m = 2 grows to a substan- 
tial amplitude with time. At about t ~ 4 Myr, it enters 
exponential growth — this explains the appearance of 
the gaseous bar at this stage. Clearly, this bar-like mode 
has formed early in the evolution, when the gas compo- 
nent does not dominate the potential at any radius, and 
when the global stability parameter a < 0.34 (see Sec- 
tion 12.11) . This mode is stimulated by the overall mass 
distribution. The exponential growth at later time hap- 
pens exactly when and where the gravitational potential 
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Fig. 9. — The temperature variation with the density in the collapsing flow of model B at four different times: Upper left: t ~ 1 Myr, 
Upper right: ~ 2 Myr, Lower left: ~ 4.61 Myr, and Lower right: ~ 4.70 Myr. Only cells at spherical radius R > 2 X 10 pc arc shown. 
The colors represent the frequency of cells in the respective mass range (right scale) . Note the appearance of the radial shock which moves 
to lower densities with time. 



of the gas becomes the dominant one. Thus the bar-like 
mode appears to be driven initially by the shape of the 
overall potential, but subsequently runs away due to the 
self-gravity of the gas. 

The temperature variations with density for model B 
are shown in Fig. [SJ At the start of disk formation, a 
radial shock forms at r ~ 2 pc and then gradually prop- 
agates out to 10 pc. At this time the central runaway is 
triggered. As shown by Fig. |H1 the posthock gas rapidly 
cools down to low T in the postshock region because of 
the high density. Within the growing disk, r <, 10 pc, we 
observe an increasing range inT^1.5xl0 3 — 10 4 K and 
p ~ 10~ 17 — 10~ 21 gcm~ 3 . The upper limit to the tem- 
perature is dropping with density. The central runaway 
at <^ 1 pc narrows the T-range to mostly T ~ 3 x 10 3 K 
but increases the range in p ~ 10~ 15 — 10~ 10 gcm~ 3 . 
The maximal postshock T in the disk is slightly increas- 
ing with time, which reflects the increasing shock-impact 
velocity of the gas which comes from progressively larger 
distances. 



4.2. Interplay between off-center fragmentation and 
central runaway 

As we noted in Section 14.11 models A - C exhibit a 
central runaway collapse, while models D and E do not 
show it — not at t = 4.7 Myr nor at anytime later on. We 
now describe the evolution of the representative model D 
and analyze the reasons for its differences from model B. 
We shall discuss additional models as well. 

The first stage of collapse in model D proceeds simi- 
larly to model B, but leads to a disk within the centrifu- 
gal barrier that has a much lower surface density at the 
center. Figure [2] shows that this disk possesses a density 
core that grows to ~ 2 — 3pc at t = 4.7 Myr — the time 
of the onset of the central runaway in model B, in sharp 
contrast with this model. In model E, the density core 
extends to ~ 8pc at this time. When the gas distribu- 
tions in models B and D are compared before the central 
runaway in model B, at t ~ 4.0 yr, the core density in 
B is higher by about 2 orders of magnitude than in D 
(Fig. 13]), and by even more compared to model E. 

Clearly, a more cuspy DM distribution leads to a more 
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Fig. 10. — (a) Face-on density-weighted projection of density of the gas disk forming from collapsing gas at t = 13.2 Myr in Run D. 
Box size is 100 pc and its vertical extent Az = 2 pc. (b) Edge-on projection of of the same disk. The disk edge is clearly visible, as well as 
the torus supported by turbulent pressure and delineated by shocks. The anisotropic density distribution is beginning to take shape and 
extends well beyond the torus. 

centrally concentrated gas distribution. This has a clear 
dynamical consequence: the dynamical timescale prior 
to the central runaway is shorter in more cuspy poten- 
tials. We observe that, at early times of the simulations, 
the m = 1 mode is the dominant mode for all the models. 
This mode is damped faster in models with more cuspy 
DM distributions, via dynamical friction, presumably of 
the gas against the DM distributions. The m = 2 mode, 
which is responsible for the central runaway collapse, can 
only grow after m = 1 has decayed significantly (Fig. [8]). 
Therefore, our models A - E form a sequence along which 
the dynamical timescale becomes longer. Any dynami- 
cal instability will have a tendency to increase its am- 
plitude more slowly along this sequence. This includes 
the growth of m = 1, 2 and 3 modes that we observe, in 
different combinations, in these models. 

Further evolution of model D shows the growth of the 
disk behind the radial shock. The central density in- 
creases by about an order of magnitude, but still falls 
short of the beginning of the central runaway in model B 
(Figure [3]). This disk growth is related to gas with an in- 
creasing angular momentum j , and steadily growing free- 
fall velocity, arriving from larger r. After ~ 6 — 7 Myr, 
the turbulent pressure in the shocked gas substantially 
exceeds the thermal pressure behind the shock. This 
leads to the formation of a geometrically-thick torus with 
a growing surface (and volume) density (Fig. [TU)) . Fig- 
ure [TT] displays the formation and evolution of this torus 
at r ~ 20 - 100 pc. The density profile at t = 13.2 Myr 
shows that a significant fraction of the collapsing gas has 
stagnated in the torus. At this time, the gas in this re- 
gion is mainly in circular motion — the radial motion 
essentially ceases and the turbulent motions decay. The 
torus becomes azimuthally inhomogeneous, and increas- 
ing density leads initially to mild off-center fragmenta- 
tion, but most of the fragments are immediately sheared 
and destroyed (e.g., Fig. H0k). Finally, at t ~ 13.4Myr, 
one of the fragments begins gravitational collapse and 
the simulation is stopped (Fig. 1X2"]) . Model E is similar 
to D, with off-center fragmentation in the torus occurring 
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Fig. 11. — Evolution of the gas density profile, averaged over 
spherical shells, after t = 4.7 Myr, in model D. Comparing with 
Figures [T] and \3\ the disk gets bigger and its central density in- 
creases. Note the formation of t he t orus outside the disk, which 
shows mild fragmentation in Fig. 1101 



at the same time, t ~ 13.4 Myr. 

We have followed the evolution of these fragments in 
models D and E. It differs substantially from the evolu- 
tion of the central runaway (i.e., collapse of the central 
fragment) in models A - C. The typical mass of the frag- 
ments is about 10 4 Mq. Each of the fragments in D and 
E collapses and exhibits a fission into two fragments. We 
did not follow their evolution further. 

The evolution of model D thus diverges from that of 
model B. The appearance of fragments perturbs the cen- 
tral region of the gas disk. We observe that these per- 
turbations lead to a loss of symmetry in the disk as well 
as displacing its center of mass, depending on the num- 
ber of fragments and their azimuthal distribution. The 
formation of an even-mode bar is suppressed. Therefore, 
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Fig. 12. — Density-weighted projection of density of (a) model D, and (b) model E, in the equatorial plane with a thickness Az = 2pc, 
show the gas disk forming from collapsing gas at t = 13.4 Myr, when the tori in models D and E exhibit fragmentation. Box size is 100 pc. 



the effect of fragmentation in the torus is that the bar 
instability is damped. 
This result reveals competition between two timescales 

— that of the central runaway and that of off-center frag- 
mentation in the models. The development of the bar 
cascade and of off-center fragmentation are facilitated by 
different instabilities. Gas inflow, of course, drives both 
instabilities, but the DM plays an important role by regu- 
lating the timescale of the central runaway, as well as the 
decay timescale for the m = 1 mode. Gaseous bar for- 
mation is a global instability in the gas disk and appears 
to be triggered by both the gas gravity and the global 
distribution of the collapsing matter. In more centrally- 
concentrated disks, this instability happens earlier. On 
the other hand, the fragmentation in the torus is deter- 
mined by the local surface density and the local pressure 

— this is a local instability. If the gas turbulent veloci- 
ties decay first in the torus (i.e., in the outer disk), the 
torus fragments first. The competition between these 
t wo instabilities, local and globa l, has been investigated 
in lBegelman fc Shlosmanl (|2009f ). 

To summarize, we find that the DM density profile, 
specifically its cuspincss, plays the role of the discrimi- 
nator between models that experience the second stage 
of gravitational collapse past the centrifugal barrier and 
models that do not exhibit this dynamical stage. This di- 
chotomy is related directly to the surface density profile 
of the growing disk within the centrifugal barrier. But, 
as we have shown in Section |4~TI it also depends on the 
asymmetry of the outer mass distribution that develops 
in the collapsing matter, triggering the bar cascade. Es- 
sentially, this corresponds to the dynamical decoupling of 
the gas from the underlying DM potential. During the 
initial stage of gas collapse, the DM potential dominates 
at all radii and the baryon density is lower than the DM 
density everywhere. However, since the DM potential of 
model B is more cuspy than that of model D (Table [1] 
and Fig. [2]), its gas contracts more, and the surface den- 
sity of the forming disk is substantially higher. The disk 
continues to grow as a result of ongoing accretion and 
the settling down of the turbulent gas. We, therefore, 
turn to the details of this turbulent dynamics. 



4.3. The role of turbulence 

ENZO has been extensively tested to handle supersonic 
turbulent motions in a uniform density background (e.g., 
iKritsuk et al.ll2007t iPadoan et al.l l 2009f ) and in stratified 
densities ()Kritsuk et al.ll2011| ). 

Figure [5] shows the radial and tangential velocity pro- 
files of the flow in Mach numbers, calculated using the 
velocity and temperature maps for model B. The right- 
hand maximum corresponds to the initial stage of the 
gravitational collapse and exceeds A4 ~ 3 — 5. The left- 
hand maximum reflects the accelerated runaway in the 
central region and exhibits the same range of Mach num- 
bers. Both the radial and tangential flows are clearly 
supersonic. 

The upper frames of Figure [13] show the velocity field 
for model B within the central 10 pc at the time of the 
central runaway. The face-on disk displays spiral shocks 
in the gas, while the velocity field is dominated by ro- 
tation. At the same time, the edge-on disk shows a 
turbulent velocity field and a number of eddies. The 
geometrically-thin disk appears substantially puffed-up 
behind the radial shock. Because the temperature maps 
reveal that the shocks are nearly isothermal, the concur- 
rent increase in the vertical thickness of the disk can come 
only from the turbulent flow behind the standing radial 
and vertical shocks. Indeed, our estimates of thermal 
pressure gradients in the z-direction reveal that thermal 
pressure gradients are insufficient to puff up the disk to 
the extent seen in the upper right frame of Figure 1131 

While turbulence is notoriously difficult to define, we 
follow the definition which relies on the vorticity w = 
V x v, and its cross product with the velocity field, i.e., 
the inertial vortex force v x V x v. To quantify the 
turbulence within the disk, we have followed the vor- 
ticity field within the computational box (lower frames 
of Fig. [T3|) . As expected, on larger spatial scales the ve- 
locity field is irrotational due to the relaxed initial con- 
ditions used in our simulations. As the inflow velocity 
grows with time and the velocity field becomes less reg- 
ular, the vorticity increases and exhibits a discontinuity 
at the standing shock which envelops the disk. The post- 
shock flow shows a sharp increase in the vorticity, which 
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Fig. 13. — Top: Face-on (left) and edge-on (right) projections of the velocity field in the collapsing flow of model B, on scales of 20 pc. 
Arrows show the direction of the flow only — the velocity values are given by the color palette. While the face-on flow is dominated by 
the rotational component (at larger radii), the edge-on slice is dominated by turbulence. Positions of radial and surface shocks are clearly 
delineated by sharply increased turbulence. Bottom: Face-on (left) and edge-on (right) slices of the vorticity field, w, in the collapsing flow 
of model B on 20 pc scales. Vorticity is generated by the oblique shocks that delineate the disk, by spiral shocks within the disk, and by 
the central runaway. 



decays toward the equatorial plane. In the postshock re- 
gion, the turbulence is transonic (Table [2]) . It provides 
support for the vertical disk structure. Its decay, when 
moving radially-inward from the shock, results in the 
sharp decrease of the disk vertical thickness at smaller 
radii. 

The vorticity in the disk appears to be driven by the 
spiral shocks and by the bar-like perturbation at the cen- 
ter (at later times) — the spiral arms around this pertur- 
bation are turbulent as well, as can be seen in the face-on 
disk region. The central region, which experiences the 
runaway collapse, exhibits the highest vorticity. 

We have also sampled the turbulent velocity field (in 
terms of the Mach number) on smaller spatial scales, at 
the time of the central runaway, using spherical sampling 
volumes whose positions in the disk plane at radii r c arc 
given in Table [2] The evolution of the turbulent velocity 



TABLE 2 
Sampling the Turbulent Velocities in model B 



Center at r c 


(pc) 


Spr 


ere Radius (pc) 


M 


7.2 X 10~ 3 






7.0 X 10~ 3 


1.85 


7.5 x 10~ 3 






7.2 x 10~ 3 


1.62 


1.0 x lO" 1 






5.0 x 10~ 2 


0.72 


1.0 






0.3 


0.48 


2.5 






1.0 


0.63 


5.0 






2.0 


0.94 



Note. — r c is the distance from the rotation axis of the center 
of the sampling sphere ( see text for additional explanations). Last 
column: the estimated Mach number within the sampling sphere. 

within the central sphere (i.e., first line in Tabled]) is also 
displayed in Figure HH 

An alternative method for measuring the properties of 
supersonic turbulence is the PDF of the gas density. As 
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Fig. 14. — Evolution of turbulent velocities in the central run- 
away region, in units of Mach number, in model B. The region 
sampled consists of a sphere positioned in the equatorial plane at 
r c = 1,500AU from the center, with a radius of 1,450 AU. This 
corresponds to the first line in Table [2] 
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Fig. 15. — The volume-averaged PDF of the gas density as a func- 
tion of logio p measured during the central runaway at t = 4.7 Myr 
for model B and sampled with ~ 1, 500 AMR cells. The sampling 
shows the PDF of the central shell of radius 20 AU-200 AU (blue 
histogram) and the central sphere of 20 AU (green histogram) . 
Shown also are the lognormal fit (black) with dispersion a ~ 1.52.M 
for the blue histogram using Equation [2] and the power law tail 
(red) for the green histogram with a slope of ~ — 1. The collapsing 
gas has been sampled at the resolution of 1 AU and the density 
fluctuations extend over 7.5 decades. The average baryon density 
in the sample spherical shell is < p >~ 4.7 X 10~ 14 gem -3 . 

discussed in Section^ the density PDF for homogeneous 
supersonic turbul ence is expected to follow a lognormal 
distribution (e.g. , IVazauez-Semad eni 199 41 lPadoanlll995t 
IPadoan fc Nordlundl 120021 : IKrumholz fc McKed I2005h . 
To estimate the PDF in different regions of model B, 
we measure the grid densities within sampling spheres 
centered at various points of interest. The sampling of 
the central region is centered at (0, 0, 0). For each sam- 



pling center, we choose the radius of a sphere to enclose 
~ 1,500 AMR cells. These spheres sample different dy- 
namical states of the collapsing gas: the central runaway, 
the outskirts of the disk, the transition region between 
the disk and the halo, and the gas at large. We carefully 
choose the size of these spheres to sample a large enough 
number of grid cells to obtain reasonable statistics, with- 
out mixing different dynamical regions in a single sphere. 

Figure [15] shows the volume-averaged density PDF in 
the central region of model B during the central runaway, 
at t = 4.7 Myr. We are unable to fit a single analytic 
lognormal PDF, for a given mean density (Equation [2]) , 
to the entire density range. Instead, we use a combina- 
tion of a lognormal and a power-law distribution, with 
the power-law tail dominating at higher densities. At 
lower densities, comparison between the measured den- 
sity PDF and the analytical lognormal PDF shows ex- 
cellent agreement, with the lognormal PDF fit extending 
over 4 decades in p. For logp > —12, the best fit is 
a power law tail with the slope of ~ — 1 for about 3.5 
decades in log p. No such tail is detected for other loca- 
tions of sampling spheres. 

Such a power-law tail has been observed previously 
in two-dimensional simulations of homogeneous, super- 
sonic hydrodynamical turbulence in its early stages, be- 
fore the format ion of sel f-gravitating clumps, i.e., before 
fragmentation (|Scalo et alJ|1998J ) . This stage is similar to 
the present models, albeit with an important difference 
- our central region is also in a state of a supersonic 
gravitational collapse and therefore has devel oped a sub- 
stanti al density gradient. In comparison, the lScalo et al.1 
(|1998| ) power- law tail has a measured slope of ~ — 2 be- 
fore the self-gravitating clumps have formed, and a slope 
of —1 when these clumps are present. In this respect, 
the central runaway in our model s agrees well with the 
similar dynamic state in the gas of iScalo et al.l (J1998J ) - 
it corresponds to the gravitational collapse of a "frag- 
ment." Such a power law is predicted for the solutions 
of Burgers equation (which is pressur eless) at high den- 
sities (e.g., iGotoh fc Kraichnanlll993T) . Power- law tails 
have also been observed in three-dimensional sim ulations 
(e.g., iFederrath et alJ[200l iKritsuk et al.j[2^Tll) . 

In a more recent work, IKritsuk et al.l (120111 ) has tar- 



geted isothermal supersonic turbulent flows in the pres- 
ence of gas self-gravity, for the purpose of determining 
the mass density PDF. A random force has been used to 
drive the turbulence on large spatial scales, in contrast 
with our models where turbulence is driven by gravita- 
tional collapse only. Despite these differences, we are in 
agreement that the power-law tail in the PDF appears 
at the time of the central runaway, when and where the 
local gravity in the gas becomes important, albeit the 
slopes are different a t the end of the simulations: —1.7 for 
IKritsuk et al.l (|2011h and —1 for our models A - C. Sam- 
pling away from the central runaway site in our models 
does not show the power-law tail, indicating that there 
are no other self-gravitating fragments within the com- 
putational box. 

Why is the power-law slope shallow e r in o ur simula- 
tions, i.e., —1.7 vs —1? IKritsuk et al.l (J2011I ) comment 
that a shallower, —1 slope does appear at high densities 
and associate this with the mass pile-up resulting from 
dynamically important angular momentum in the region. 
Recall that the central runaway in our simulations is an- 
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gular momentum-dominated, as we show in Section 14.11 
We have tested the origin of this power-law PDF by 
sampling the region with concentric spheres having pro- 
gressively smaller radii. We find that the range fit by 
the lognormal PDF shrinks along this sequence, primar- 
ily because of cutting off the lower densities, while the 
high-density end remains untouched (compare the green 
and blue histograms in Figure [TS]). The high-density 
end of the power-law PDF, therefore, corresponds to the 
very high central density in our simulations, and the den- 
sity stratification in the central region is responsible for 
the —1 slope power-law PDF. Clearly, understanding of 
the power tail in the PDF requires additional theoretical 
analysis. 

In our models discussed here, the initial conditions are 
simplified and do not include turbulent motions. This de- 
lays the onset of turbulence, which takes time to fully de- 
velop in the outer part of the flow. In more realistic mod- 
els, the infalling gas is expected to be already partly tur- 
bulent. However, o ne can also argue i n favor of initially 
lamin ar flow (e.g., lAbel et all 120021 : iBromm fc Larsonl 
120041) . As the infall develops, turbulence sets in spon- 
taneously and greatly increases after the gas has crossed 
the standing shock enveloping the disk within the cen- 
tral few parsecs. Turbulence which has developed during 
gravitational collapse has been noticed by other authors 
(iLevine et al.ll200l IwTse et al.|[200l iRegan fc Haehneltl 
I2009T) . and its impact on t he gravitational stability of th e 
flow has been analyzed bv lBegelman fc ShlosmanT(|2009D . 
As discussed in Section [5J the effect of turbulence de- 
veloping during gravitational collapse is to suppress gas 
fragmentation. 

The absence of fragmentation in our simulations of col- 
lapsing flow is evident in Figure [1] which provides a snap- 
shot of the face-on disk in model B. A simple estimate 
based on the sound speed of the gas within a 2 x 10 8 M 
DM halo shows that the ratio of the gas mass within 
some spherical radius R within this halo to the Jeans 
mass at this temperature and density is of order ~ 3 — 4. 
Because the Jeans mass depends on the rms velocity to 
the third power, this ratio will decrease to ^ 1 in a flow 
with transonic turbulent velocities. 

4.4. Randomizing gas motions: model Dmod 

We have rerun model D with less organized baryonic 
initial momenta (Section 14.21 and Table [I]). The purpose 
of this run is to introduce more realistic initial conditions 
expected in the cosmological context. In model Dmod, 
the inner spherical shells have their j slightly decreased 
and their orientation randomized. To keep the total an- 
gular momentum unchanged, the outer shells have to 
compensate for this and their j has been slightly in- 
creased. The large-scale evolution differs little from that 
of model D. Namely, the centrifugal barrier has moved 
inward somewhat. The disk forms with higher surface 
density and is geometrically thicker. The disk thickness 
appears independent of r, unlike in Figs. [TU1 and [T51 We 
also observe that the position of the disk's equator is 
less stable relative to the equatorial plane of the DM 
halo, moving periodically along the z-axis, and that the 
density peak in the disk experiences some weak, low- 
amplitude m = 1 perturbations. 

Probably the most important difference is the dramatic 
weakening of the outer torus. This evolution has a dual 



effect — the timescale for the onset of the central run- 
away is shortened to 4.5 Myr, and the off-center fragmen- 
tation does not materialize because the surface density 
of the torus is very low. As a result, unlike in model D, 
model Dmod experiences a 'classical' central runaway, 
similar to models A - C. 

5. DISCUSSION AND SUMMARY 

We have studied some of the physical processes that 
can operate during the first stages of SMBH formation 
at high z within the direct collapse paradigm. The ini- 
tial conditions used in this work have been substantially 
simplified and consist of an isolated, responsive, spher- 
ical DM halo of M v ; r ~ 2 x 10 8 M Q and a virial radius 
of ~ 1 kpc, having a spin parameter A ~ 0.05; embedded 
baryons with an average cosmological fraction; a uni- 
versal angular momentum distribution; and nonsingular 
isothermal density profiles for DM and gas. We limit 
ourselves to the primordial composition and the absence 
of molecular cooling. 

In a set of high-resolution numerical models we have 
only varied the size of the DM density core, i.e., the re- 
gion where the DM density is constant. The collapsing 
flow is followed down to spatial scales ~ 10 -4 pc (20 AU), 
over a dynamic range of ~ 7 decades. In these simula- 
tions, we have assumed that the inflow is optically thin. 
This is tested a posteriori — the flow remains optically 
thin for electron scattering and free-free opacities. The 
Lya photons produced in the inflow have a large optical 
depth, but our estimates show that they do not have an 
effect on the dy namics of the gravitati onal collapse, in 
agreement with IRees fc Ostrikerl (|1977t ). Our modeling 
has been stopped at the approximate radius where the 
physical conditions of radiation transfer are expected to 
modify the flow substantially. The inner boundary of the 
flow should also be investigated in connection with ad- 
ditional physical processes, e.g., magnetic torques. We 
expect the weak magnetic field advectcd across the virial 
radius to be amplified by the dynamo, and the compres- 
sion to become significant here. MHD processes may also 
trigger an outflow — all this remains outside the scope 
of this work. 

Two processes can dramatically affect the outcome of 
direct collapse and prevent a seed SMBH with substan- 
tial mass from forming: efficient fragmentation and an 
angular momentum barrier. The former process can lead 
to star formation, which will sap the available mass sup- 
ply and can also expel baryons from the DM halo al- 
together, e.g., by supernovae feedback and massive stel- 
lar winds. The latter process can stop the collapse at 
the centrifugal barr ier, which has b een estimated to lie 
at ~ 0.1i?vir (e.g.. IMo et al.l I1998D . Our results show 
that efficient fragmentation has been damped by the 
development of supers o nic tu rbulence, as suggested by 
IBegelman fc Shlosm an (2009). On the other hand, we 
find that the angular momentum is not conserved within 
the centrifugal barrier — both the outer baryonic col- 
lapse and increasing self-gravity in the interior flow trig- 
ger the growth of the m = 2 bar-like mode, which chan- 
nels the gas inward. With more realistic initial condi- 
tions, the typical triaxiality of the DM halo will amplify 
the non-conservation of angular momentum. 

Overall, we find that direct collapse within DM halos 
depends on the competition between two timescales - 
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that of the central runaway within the centrifugal bar- 
rier, and that of off-center fragmentation. If we take 
a broader view, the centrifugal barrier appears to be a 
typical rather than exceptional feature of such a collapse. 
However, in all models it is situated initially much deeper 
than anticipated, at ~ 1 pc compared to the expected 
~ 100 pc for such DM halos. In models with larger DM 
cores, the central runaway time is delayed, and the ra- 
dial shock has time to propagate much farther out to 
~30-40pc. 

Why does the centrifugal barrier lie so much deeper 
than anticipated? The reason for this is the inside-out 
development of the collapse, where only the inner gas 
has time to reach the barrier. For example, in model 
B, the centrifugal barrier and the radial shock, which 
delineate it, form at about 1 pc and advance to about 
10 pc by the end of the simulation. We note an impor- 
tant detail — our simulations extend over a timescalc 
which is much shorter than the global free-fall time, 
- (3tt/32G < p >) x / 2 - 80Myr, in these DM halos. The 
central runaway is triggered within the first 2% — 18% of 
this time, and lasts for < 1 Myr. This explains why the 
supersonic turbulence did not decay in our models and 
why the fragmentation process is so inefficient. 

The most intriguing characteristics of the central run- 
away are that it is self-similar and disky. This means that 
the angular momentum is dynamically important down 
to the optically-thick boundary of about 10 AU. The col- 
lapsing region is partially supported by rotation in the 
radial direction and pressure supported in t he vertical 
direc tion. Published models of SMS, (e.g., iBegelmanl 
l2010f ). do assume a degree of rotation in order to en- 
sure stability, but it is well below the dynamically im- 
portant rotation encountered near the inner boundary of 
our models (e.g.. iMontero et ail 12012ft . It is, therefore, 
natural to assume that the transition from a rotationally- 
supported entity to a pressure-supported SMS happens 
close to this boundary. However, the details of this tran- 
sition are completely unclear. Moreover, a possibility 
exists that the runaway collapse dominated by j will by- 
pass the state of a hydrostatic equilibrium, that ther- 
monuclear fusion will not play a major role, and that 
the collapse will proceed directly to forming the SMBH 
horizon. 

To properly follow up the gravitational collapse it is 
critical to resolve the centrifugal barrier and the associ- 
ated radial shock, i.e., to have a spatial resolution of a 
fraction of a parsec. At lower resolution the evolution 
can diverge from the one we observe. 

Usage of a randomized j (R) orientation leads to the fol- 
lowing corollaries. The centrifugal barrier moves slightly 
inward, the disk becomes somewhat smaller -- these 
changes do not appear significant. But the disk surface 
density increases substantially. As a result, the central 
runaway in model Dmod happens after ~ 4.5 Myr. The 
corresponding model D exhibits off-center fragmentation 
instead. The largest difference between these models is 
the absence of the massive torus that dominates the outer 
disk in D. Only a trace of this configuration remains in 
Dmod, and it is stable against fragmentation. So more 
random initial conditions move model D into the 'main- 
stream' of models A - C. Clearly, cosmological initial 
conditions will show the most realistic solution. 



Another important requirement is to resolve the super- 
sonic turbulence which develops in various parts of the 
collapsing flow, and especially behind the radial shock 
and during the central runaway. The relative absence 
of turbulence between the shock and the virial radius 
comes from the quiescent flow in this region — a direct 
consequence of our initial conditions of an isolated DM 
halo. In more realistic cosmological initial conditions the 
laminar flow around i? V ir may be already turbulent. 

We have analyzed the density PDF in the central run- 
away and found that it is not the lognormal PDF typ- 
ically encountered in simulations of non-self-gravitating 
isothermal supersonic flows. The PDFs in models A - 
C consist of the usual lognormal part as well as a high- 
density power-law part. The slope of the power-law is 
found to be ~ —1 at the time and position of the central 
runaway. Limiting the sampling of the density fluctu- 
ations to a smaller region close to the very center, we 
find that the lognormal PDF fades away but the power 
law part remains intact (Fig. II 5|) . The lognormal density 
PDF extends over 4 decades in density and the power- 
law extends over 3.5 additional decades at higher density. 
Comparison with two-dimensi onal hydrodynam ical sim- 
ulations with gas self-gravity (jScalo et al.lll998l see also 
Kritsuk et al. 2011 for 3-D simulations) confirms that 
the formation of self-gravitating clumps in the presence 
of supersonic turbulence depends on the position of the 
velocity sampling and shows the same structure of a log- 
normal + power-law PDF. 

5.1. Estimating the seed SMBH mass range 

We now turn to estimating the seed SMBH masses. Ta- 
ble [3] shows the onset time of the central runaway, t co u , 
in models A - C, increasing from 2.1 Myr (A), to 4.7Myr 
(B), to 8.7Myr (C). Models D and E do not show central 
collapse before we observe off-center fragmentation. We 
have calculated the radial dependence of the mass accre- 
tion rate, M(R), for all models (e.g., Fig.[6]for model B). 
In all cases, the central runaway extends radially over a 
fraction of the central disk. We use the radial mass accre- 
tion rate profiles to estimate the baryon mass that par- 
ticipates in the central runaway. In Figure [6] for model 
B, as well as for models A and C, we observe that these 
baryons are the ones located within about the half-radius 
of the disk at the runaway time, as given in Table Note 
that the disk radius is given by the position of the radial 
shock at the centrifugal barrier. This also corresponds 
to the global minimum of M(R) at that time — baryons 
within this radius effectively decouple from the DM back- 
ground and are dumped onto the center. Baryon masses 
which are part of this breakaway are also given in Table [3] 
and range between M co ii ~2x 10 4 M Q and 6 x 10 5 M Q . 

We now attempt to assess the validity of these esti- 
mates. Plotting M co ii as function of the onset of the cen- 
tral runaway collapse time, £ C oii, gives a nearly perfect 
log-linear dependence, logM co u ~ t C oii, namely, 



log (M coU /M Q ) = a(*coii/l Myr) + b, 



(3) 



where a ~ 0.18 and b ~ 3.95. On the other hand, £ co ii 
depends linearly on the size of the DM density core in 
models A - C, i?dm, given in Table [TJ Assuming that 
M co u has a direct relationship to the mass of the SMBH 
seed, M», models with larger £ co n, which lead to larger 



Runaway Collapse 



17 



TABLE 3 
Parameters of the Central and Off-Center Runaways 



Models 


icon (Myr) 


-Rcoll (pc) 


iWcoll (Me) 


A 


2.1 


2 


2 X 10 4 


B 


1.7 


5 


8 x 10 4 


C 


8.7 


10 


6 x 10 5 


D 


13.4 


off-center 


- 


E 


13.4 


off-center 


- 


Dmod 


4.5 


5 


2 x 10 5 



Note. — i co n - onset time of the central runaway or of the off- 
center fragmentation; R co ii — initial radius of the central runaway; 
-Mcoll - baryon mass participating in the central runaway. 



M co \\, should also result in larger M.. However, an up- 
per limit on i co n appears to come from the condition for 
off-center fragmentation in the torus — this limit comes 
from models D and E, which both show off-center frag- 
mentation at ~ 13.4 Myr. Hence £ co ii Ss 13.4 Myr, which 
is rather a conservative estimate as the fragments will 
need some time to affect the central runaway dynamics. 

The upper limit of 13.4 Myr intersects the £ co ii(-Rs) line 
at R s ~ 1.22 pc. Models with larger R^ m should exhibit 
off-center fragmentation. We test this on models D and 
E — both lie to the right of 1.22 pc (Table [TJ. So our 
simplistic argument has passed the first test successfully. 
What have we learned from this reasoning? 

The central runaway drains baryons within the radius 
~ -Rcoii, when the gas accumulation inside this radius 
roughly exceeds that of the DM. The collapse time can 
be estimated roughly as ~ 2 — 3 x iff, where tg is the 
local, i.e., within ~ i? C oib free-fall timescale. Thus the 



collapse timescale is ~ 3 x W 6 R coll 10 ^ ( 



Af_„n a yrs, where 



#0011,10 = ftcoii/10pc and M co u i6 = 7\/ co n/lO b M . One 
should consider that baryons inside i? C oii can be replen- 
ished, in principle, as the material flows in across the 
radial and surface shocks. This would determine a char- 
acteristic timescale which may be an order of magnitude 
above the estimated collapse time. 

Probably the most intriguing consequence of this argu- 
ment is the emerging mass range for the SMBH seeds. If 
a large fraction of the overall inflow goes into formation 
of the SMBH seed, we can extrapolate log M co \\ ~ t co n to 
obtain the maximal M co \\ ~2x 10 6 M , which is about 
10% of the amount of baryons in DM halos of interest, 
M vir ~ 1 - 2 x 1O 8 M . So the mass range for SMBH 
seeds appears to be 2 x 10 4 M Q <, M. <, 2 x 10 6 M Q . If, 
in addition, the size of the flat DM density core corre- 
lates with the halo virial radius, the mass of the SMBH 
seed is expected to correlate with the DM halo mass, at 
the time of formation. This also hints at the possible 
correlation between the DM halo mass function and the 
SMBH seed mass function. 

Our models relate the properties of SMBHs formed 
through direct collapse to the sizes of the flat density 
cores of DM halos. Pure DM simulations (e.g., NFW) 
have claimed universal density profiles with a cusp, while 
observations hi nt rather at the existen c e of flat den- 
sity cores (e.g ., iFlores fc Primackl 119941: Ide Blokl 120051: 
IPrimackll2009l for review). A possible explanation for the 
flattening of NFW density cusps, appea l ing to the action 
of clumpy baryons (lEl-Zant et al.ll200lL I2004D. has been 
verifi ed in numerical simulations (|Romano-Diaz et al.l 
|2008[) . Other solution within the CDM paradigm rely on 



baryon energy feedback (e.g.. lMashchenko et al J 120061) . 

The size of the DM density cusp in the NFW profiles, 
and, therefore, the size of the DM density core replacing 
the cusp, strictly correlate with the halo virial radius. 
This assumption is probably overly optimistic, and relies 
heavily on the fragility of the cusp due to its thermo- 
dynamic improbability (|E1-Zant et al.|[200l| ). Neverthe- 
less, we point out that such a correlation will lead to an 
SMBH mass which initially depends linearly on the DM 
halo mass. In this case the SMBH seed mass can be in- 
ferred from TableGHto lie at M. ~ 10 6 M Q for a DM halo 
of M v i r ~2x 10 8 Mq and i? v i r ~ 1 kpc, which will have 
a DM density core of i?dm ~ 1 pc. This is close to the 
upper limit on M m we have estimated above. If, however, 
i?dm does not correlate with i? v ir, the above estimate of 
an SMBH mass range 2 x 10 4 M <, M. <, 2 x 10 6 M Q 
appears to be more realistic. 

The above conclusions will be modified if a substantial 
amount of the collapsing baryons is expelled via some 
feedback or wind mechanism, effectively decreasing the 
peak accretion below its nominal value of M ~ 1 — 2 M . 
moreover, the proposed range of M m is attributed only 
to a single cycle in the accretion process, by which we 
mean one central runaway resulting in the formation of 
the SMBH seed. The conditions leading to the second 
cycle will differ because of the existence of the SMBH. 
However, it is not clear if the mechanical and radiative 
feedback from this seed will have an effect on the next 
runaway, i.e., on the 2nd cycle. One can envision that 
the feedback is directed along the rotation axis, while the 
next central runaway proceeds in the equatorial plane. 

Finally, we note that while our initial calculations of 
SMBH formation in the direct collapse scenario have em- 
phasized some interesting outcomes, a long list of issues 
to be resolved remains. One such issue is whether molec- 
ular hydrogen can affect the outcome of this process by 
inducing fragmentation in the collapsing gas. Nearby 
stars can contribute to the UV backgrou nd which have 
an ad verse effect on H2 formation (e.g., Diikstra et al. 
l2008f ). In this work we assume that the UV back- 
ground will damp H2 formation and, therefore, will 
maintain the gas temperature not much below T gas ~ 
10 4 K. Several physical mechanisms have been proposed 
that can support th i s state of the gas (lOmukail I200H 
Oh fc HaimarJ [20021 IsiTaans fc Silld I200fit IShang eTall 



2010L ISchleicher et all I2010t lLatif et all I2011J L Imple- 
mentation of radiative transfer calculations to study the 
details of this process is a nex t logical step. We note that 
iBegelman fc Shlosmanl ([2009D have argued that fragmen- 
tation will be suppressed even if the cooling floor moves 
substantially below 10 4 K. This happens because the flow 
becomes much more supersonic and the fraction of frag- 
menting gas decreases with increasing M. . Alternatively, 
if the collapse happens inside more massive halos, the ra- 
tio of the virial velocity to the sound speed will increase 
and lead to the same result. 

Although our initial conditions have been motivated 
by the current cosmology framework, they are signifi- 
cantly simplified and idealized. Owing to this simplifi- 
cation, the simulation results can qualitatively demon- 
strate the physical processes that work but cannot 
quantitatively predict the physical timescales discussed 
here. Varying the DM halo profile, gas density pro- 
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files, and angular momentum distribution can affect the 
bar formation timescale and the torus fragment ation 
timescale. For example. iKoushiappas et al] (|2004l ) and 
lLodato fc Nataraianl (J2006I ) suggested that early SMBHs 
tend to be formed in halos with a low angular momen- 
tum. It will be interesting to predict the environments of 
early SMBHs formed through direct collapse, using full 
cosmological simulations with many different halo condi- 
tions. 
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